if (class(model) %in% "character") {
  model <- ModelTBBCGEngland::read_libbi(file.path(model_dir, "libbi", "posterior"))
}

priors <- readRDS(file.path(model_dir, "data", "prior-params.rds"))
model
## Wrapper around LibBi
## ======================
## Model:  Baseline 
## Run time:  2070.558  seconds
## Number of samples:  1000 
## State trajectories recorded:  E H L P S T_E T_P YearlyAgeCases YearlyCases YearlyEAgeCases YearlyECases YearlyEPulCases YearlyEPulDeaths YearlyNonUKborn YearlyPAgeCases YearlyPCases YearlyPulCases YearlyPulDeaths 
## Observation trajectories recorded:  YearlyAgeInc YearlyHistPInc YearlyInc 
## Parameters recorded:  alpha_t_decay alpha_t_init c_eff c_hist chi_init delta epsilon_h_0_4 epsilon_h_15_89 epsilon_h_5_14 epsilon_l_0_4 epsilon_l_15_89 epsilon_l_5_14 HistMeasError kappa_0_4 kappa_15_89 kappa_5_14 M mu_e_0_14 mu_e_15_59 mu_e_60_89 mu_p_0_14 mu_p_15_59 mu_p_60_89 nu_e_0_14 nu_e_15_89 nu_p_0_14 nu_p_15_89 phi_0_14 phi_15_59 phi_60_89 rho_0_14 rho_15_59 rho_60_89 Upsilon_0_14 Upsilon_15_59 Upsilon_60_89 zeta_e_0_14 zeta_e_15_59 zeta_e_60_89 zeta_p_0_14 zeta_p_15_59 zeta_p_60_89
traces <- coda::mcmc(rbi::get_traces(model))

coda::rejectionRate(traces)
##               M           c_eff          c_hist           delta 
##        0.992993        0.992993        0.992993        0.992993 
##   epsilon_h_0_4  epsilon_h_5_14 epsilon_h_15_89       kappa_0_4 
##        0.992993        0.992993        0.992993        0.992993 
##      kappa_5_14     kappa_15_89   epsilon_l_0_4  epsilon_l_5_14 
##        0.992993        0.992993        0.992993        0.992993 
## epsilon_l_15_89        phi_0_14       phi_15_59       phi_60_89 
##        0.992993        0.992993        0.992993        0.992993 
##    Upsilon_0_14   Upsilon_15_59   Upsilon_60_89        rho_0_14 
##        0.992993        0.992993        0.992993        0.992993 
##       rho_15_59       rho_60_89       nu_p_0_14      nu_p_15_89 
##        0.992993        0.992993        0.992993        0.992993 
##       nu_e_0_14      nu_e_15_89     zeta_p_0_14    zeta_p_15_59 
##        0.992993        0.992993        0.992993        0.992993 
##    zeta_p_60_89     zeta_e_0_14    zeta_e_15_59    zeta_e_60_89 
##        0.992993        0.992993        0.992993        0.992993 
##       mu_p_0_14      mu_p_15_59      mu_p_60_89       mu_e_0_14 
##        0.992993        0.992993        0.992993        0.992993 
##      mu_e_15_59      mu_e_60_89        chi_init    alpha_t_init 
##        0.992993        0.992993        0.992993        0.000000 
##   alpha_t_decay   HistMeasError 
##        0.992993        0.992993
coda::autocorr.plot(traces)

- Evaluate Posterior Traces

plot_param(model, scales = "free", plot_type = "trace", burn_in = 0) + theme(legend.position = "none")

plot_param(model, prior_params = priors, scales = "free")

param_sum <- plot_param(model, prior_params = priors, plot_data = FALSE) %>% 
  group_by(Distribution, parameter, length) %>% 
  summarise(median = median(value), 
            lll = quantile(value, 0.025), 
            hhh = quantile(value, 0.975)) %>% 
  mutate(value = pretty_ci(median, lll, hhh, sep = ", ")) %>%
  group_by(Distribution, parameter) %>% 
  mutate(Parameter = case_when(max(length) > 1 ~ paste(parameter, 1:n(), sep = "-"),
                               TRUE ~ as.character(parameter))
         ) %>% 
  ungroup %>% 
  select(Distribution, Parameter, value) %>% 
  spread(key = "Distribution", value = "value") %>% 
  select(Parameter, Prior, Posterior)

saveRDS(param_sum, file.path(model_dir, "data", "sum_prior_posterior.rds"))
knitr::kable(param_sum)
Parameter Prior Posterior
alpha_t_decay -0.13 (-0.23, -0.04) -0.07 (-0.13, -0.04)
alpha_t_init 0.84 (0.76, 0.90) 0.79 (0.13, 0.87)
c_eff 2.42 (0.13, 4.87) 2.21 (0.16, 4.44)
c_hist 12.57 (10.14, 14.82) 12.99 (10.71, 13.16)
chi_init 0.19 (0.08, 0.30) 0.14 (0.08, 0.25)
delta 0.78 (0.70, 0.86) 0.78 (0.73, 0.80)
epsilon_h_0_4 1.53 (0.08, 4.67) 0.49 (0.20, 5.01)
epsilon_h_15_89 23.87 (1.15, 74.94) 20.29 (19.09, 57.16)
epsilon_h_5_14 3.66 (0.23, 11.82) 3.50 (0.37, 14.25)
epsilon_l_0_4 581.14 (39.24, 1657.56) 523.98 (191.49, 802.66)
epsilon_l_15_89 1106.77 (53.74, 3227.33) 98.68 (78.38, 101.18)
epsilon_l_5_14 522.01 (29.10, 1523.41) 311.46 (149.90, 1227.31)
HistMeasError 1.00 (0.63, 1.42) 1.16 (0.89, 1.36)
kappa_0_4 0.80 (0.04, 2.49) 1.13 (0.25, 1.85)
kappa_15_89 1.15 (0.07, 3.44) 0.90 (0.12, 1.98)
kappa_5_14 0.98 (0.04, 3.09) 0.34 (0.06, 1.25)
M 0.25 (0.01, 0.48) 0.27 (0.04, 0.41)
mu_e_0_14 276.37 (206.60, 336.17) 242.32 (177.20, 268.74)
mu_e_15_59 188.30 (76.36, 317.88) 237.15 (129.68, 311.23)
mu_e_60_89 35.34 (1.87, 102.13) 2.22 (0.44, 10.87)
mu_p_0_14 244.82 (156.68, 328.74) 242.42 (202.11, 261.82)
mu_p_15_59 85.76 (4.20, 250.49) 100.64 (51.90, 139.53)
mu_p_60_89 51.88 (2.09, 165.39) 72.92 (11.51, 103.42)
nu_e_0_14 0.17 (0.00, 1.40) 0.00 (0.00, 1.54)
nu_e_15_89 0.32 (0.01, 1.83) 0.33 (0.06, 1.15)
nu_p_0_14 0.11 (0.00, 0.68) 0.06 (0.05, 0.19)
nu_p_15_89 0.25 (0.02, 1.10) 0.07 (0.02, 1.10)
phi_0_14 0.58 (0.28, 1.01) 0.57 (0.49, 0.75)
phi_15_59 0.62 (0.27, 1.19) 0.42 (0.41, 0.80)
phi_60_89 0.61 (0.28, 1.09) 0.79 (0.48, 0.86)
rho_0_14 0.30 (0.27, 0.34) 0.30 (0.27, 0.32)
rho_15_59 0.65 (0.64, 0.66) 0.66 (0.65, 0.66)
rho_60_89 0.54 (0.52, 0.55) 0.54 (0.53, 0.55)
Upsilon_0_14 0.63 (0.61, 0.65) 0.64 (0.61, 0.65)
Upsilon_15_59 0.71 (0.70, 0.71) 0.70 (0.70, 0.71)
Upsilon_60_89 0.75 (0.74, 0.76) 0.75 (0.74, 0.76)
zeta_e_0_14 123.90 (55.57, 191.79) 135.23 (95.99, 154.38)
zeta_e_15_59 60.60 (3.08, 184.34) 31.35 (19.77, 77.57)
zeta_e_60_89 133.27 (61.54, 208.20) 102.43 (27.49, 153.23)
zeta_p_0_14 94.05 (16.33, 176.11) 63.57 (42.05, 116.35)
zeta_p_15_59 82.12 (2.64, 255.22) 10.10 (4.90, 206.49)
zeta_p_60_89 124.89 (20.35, 252.94) 56.52 (38.38, 308.59) 
plot_state(model, summarise = TRUE)

plot_state(model, summarise = TRUE, start_time = 59)

plot_state(model, states = c("YearlyHistPInc", "YearlyInc"), summarise = TRUE, start_time = 49) -> p2

p2

pred_obs <- plot_state(model, states = c("YearlyHistPInc", "YearlyInc", "YearlyAgeInc"), summarise = FALSE, start_time = 59, plot_data = FALSE)

obs <- ModelTBBCGEngland::setup_model_obs() %>% 
  bind_rows(.id = "state") %>% 
  mutate(time = time + 1931)

pred_obs <- pred_obs %>% 
  dplyr::filter(Average == "median") %>% 
  mutate(pred_value = pretty_ci(Count, lll, hhh, digits = 0)) %>% 
  left_join(obs, by = c("state" = "state", "time" = "time", "age" = "age")) %>% 
  select(Observation = state, Age = age, Year = time,  `Observed Incidence` = value, `Predicted Incidence` = pred_value) %>% 
  mutate(Year = Year)


sum_obs_table <- pred_obs %>% 
  dplyr::filter(Observation %in% c("YearlyInc", "YearlyHistPInc")) %>% 
  drop_na(`Observed Incidence`) %>% 
  select(-Age) %>% 
  mutate(Observation = case_when(Observation == "YearlyInc" ~ "All TB cases",
                                 Observation == "YearlyHistPInc" ~ "Pulmonary TB cases"))


saveRDS(sum_obs_table, file.path(model_dir, "data", "sum_obs_table.rds"))

age_cases_table <- pred_obs %>% 
  dplyr::filter(Observation %in% c("YearlyAgeInc")) %>% 
  drop_na(`Observed Incidence`) %>% 
  select(-Observation) %>% 
  group_by(Year) %>% 
  mutate(Age = c(paste0(seq(0, 65, 5), "-", seq(4, 69, 5)), "70-89")) %>% 
  ungroup


saveRDS(age_cases_table, file.path(model_dir, "data", "age_cases_table.rds"))
knitr::kable(sum_obs_table)
Observation Year Observed Incidence Predicted Incidence
Pulmonary TB cases 1990 3618 3448 (2734 to 3752)
Pulmonary TB cases 1991 3596 3310 (2855 to 3585)
Pulmonary TB cases 1992 3816 3297 (2816 to 3623)
Pulmonary TB cases 1993 3961 3292 (2645 to 3569)
Pulmonary TB cases 1994 3694 3167 (2690 to 3547)
Pulmonary TB cases 1995 3711 3180 (2733 to 3431)
Pulmonary TB cases 1996 3690 3048 (2616 to 3400)
Pulmonary TB cases 1997 3947 3106 (2626 to 3449)
Pulmonary TB cases 1998 3926 3091 (2527 to 3391)
Pulmonary TB cases 1999 3827 2975 (2496 to 3211)
All TB cases 2000 1803 2152 (2043 to 2521)
All TB cases 2001 1866 1949 (1948 to 2503)
All TB cases 2002 1833 2012 (1849 to 2403)
All TB cases 2003 1685 1919 (1919 to 2297)
All TB cases 2004 1776 1795 (1747 to 2201) 
knitr::kable(age_cases_table)
Age Year Observed Incidence Predicted Incidence
0-4 2000 75 2 (0 to 3)
5-9 2000 59 1 (0 to 8)
10-14 2000 75 0 (0 to 13)
15-19 2000 112 7 (3 to 12)
20-24 2000 133 15 (15 to 21)
25-29 2000 116 45 (32 to 49)
30-34 2000 147 81 (64 to 86)
35-39 2000 99 129 (108 to 163)
40-44 2000 83 175 (146 to 228)
45-49 2000 105 227 (214 to 265)
50-54 2000 116 323 (253 to 329)
55-59 2000 112 318 (318 to 382)
60-64 2000 101 326 (309 to 423)
65-69 2000 99 300 (297 to 397)
70-89 2000 371 229 (187 to 251)
0-4 2001 113 3 (1 to 5)
5-9 2001 65 1 (0 to 10)
10-14 2001 51 1 (0 to 10)
15-19 2001 130 3 (2 to 11)
20-24 2001 145 14 (9 to 17)
25-29 2001 127 32 (27 to 44)
30-34 2001 122 76 (54 to 77)
35-39 2001 128 98 (84 to 122)
40-44 2001 107 180 (153 to 198)
45-49 2001 101 215 (209 to 266)
50-54 2001 96 278 (246 to 356)
55-59 2001 91 335 (310 to 363)
60-64 2001 94 336 (306 to 379)
65-69 2001 105 307 (300 to 363)
70-89 2001 391 214 (184 to 262)
0-4 2002 102 4 (1 to 4)
5-9 2002 62 1 (1 to 13)
10-14 2002 64 1 (0 to 19)
15-19 2002 129 7 (3 to 7)
20-24 2002 148 12 (10 to 24)
25-29 2002 120 36 (27 to 47)
30-34 2002 112 59 (54 to 69)
35-39 2002 127 112 (109 to 155)
40-44 2002 98 147 (141 to 185)
45-49 2002 89 220 (200 to 254)
50-54 2002 104 282 (235 to 344)
55-59 2002 116 328 (278 to 381)
60-64 2002 88 309 (286 to 389)
65-69 2002 90 328 (273 to 390)
70-89 2002 384 214 (199 to 262)
0-4 2003 74 2 (2 to 5)
5-9 2003 46 1 (0 to 19)
10-14 2003 59 1 (1 to 13)
15-19 2003 102 6 (4 to 9)
20-24 2003 116 13 (6 to 17)
25-29 2003 137 26 (25 to 39)
30-34 2003 118 54 (54 to 73)
35-39 2003 113 98 (82 to 120)
40-44 2003 109 152 (128 to 180)
45-49 2003 83 201 (181 to 267)
50-54 2003 102 252 (246 to 305)
55-59 2003 92 282 (276 to 364)
60-64 2003 98 288 (281 to 357)
65-69 2003 94 290 (279 to 348)
70-89 2003 342 242 (199 to 298)
0-4 2004 112 4 (0 to 4)
5-9 2004 71 0 (0 to 17)
10-14 2004 81 1 (0 to 10)
15-19 2004 111 5 (2 to 6)
20-24 2004 120 14 (13 to 16)
25-29 2004 136 28 (17 to 32)
30-34 2004 129 53 (35 to 62)
35-39 2004 103 101 (67 to 126)
40-44 2004 128 132 (131 to 143)
45-49 2004 94 188 (187 to 242)
50-54 2004 103 240 (223 to 280)
55-59 2004 98 282 (281 to 342)
60-64 2004 86 297 (277 to 385)
65-69 2004 92 307 (279 to 358)
70-89 2004 312 212 (211 to 280) 
plot_state(model, states = c("YearlyAgeInc"), start_time = 59) + theme(legend.position = "none")